#----------------------------------------------------------#
#  Democratic Resilience (Aurel Croissant and Lars Lott)   #
#----------------------------------------------------------#

# Title: "Democratic Resilience in the Twenty-First Century. Search for an Analytical Framework and Explorative Analysis" #
# Authors: "Croissant, Aurel and Lott, Lars", Heidelberg University and FAU Erlangen-Nürnberg
# date: 2025-05-07
# journal: Political Studies
# DOI: 10.1177/00323217251345779
# written under "R version 4.5.0 (2025-04-11 ucrt)"

#### Preliminaries ####

R.version$version.string

# clear workspace
rm(list=ls())

# load packages
library(here)
library(tidyverse)
library(magrittr)
library(reshape2)
library(ggpubr)
library(gapminder)
library(countrycode)
library(ggrepel)
library(ggthemes)
library(lubridate)
library(dotwhisker)
library(estimatr)
library(texreg)
library(readxl)
library(haven)
library(scales)
library(zoo)
library(MCMCpack)
library(ExPanDaR)

## Install local v-dem tools package
devtools::install_local("vutils", force = TRUE)
library(vutils)

set.seed(543)

#---------------------------------------------------#
#--------------------load data ---------------------#
#---------------------------------------------------#

# set you working directory here #

## load V-Dem data V14##

vdem_full_v14 <- readRDS(file.path("data/vdem14", "V-Dem-CY-Full+Others-v14.rds"))

## load V-Party Data V2 ##

v_party_v2 <- readRDS(file.path("data", "v_party" ,"V-Dem-CPD-Party-V2.rds"))

## Mood and Satisfaction data ##

democratic_satisfaction_data <- read.csv("data/claassen/satis_est_v2.csv")
democratic_mood_data <- read.csv("data/claassen/mood_est_v5.csv")

## load political trust data from BJPS study (2025) ##

civil_mood_data <- read.csv("data/political trust BJPS study/civil_mood_est.csv")
gov_mood_data <- read.csv("data/political trust BJPS study/gov_mood_est.csv")
leg_mood_data <- read.csv("data/political trust BJPS study/leg_mood_est.csv")
parl_mood_data <- read.csv("data/political trust BJPS study/parl_mood_est.csv")
police_mood_data <- read.csv("data/political trust BJPS study/police_mood_est.csv")
polpar_mood_data <- read.csv("data/political trust BJPS study/polpar_mood_est.csv")

## store country_unit for BFA ##

country_unit <- vdem_full_v14 %>%
  dplyr::select(country_text_id, country_id, year, project, historical_date) %>%
  dplyr::filter(year >=2000)

saveRDS(country_unit, "data/bfa/country_unit.rds")


##################################################
#### Build  overall dataset ######################
##################################################

#### Create basic dataset ####

ds.basic <- countrycode::codelist_panel %>% 
  filter(year >= 1946) %>% # reduce to post-1946 period
  rename(country = country.name.en) %>% # reduce to post-1946 period
  dplyr:: select(country, year, continent, cown, vdem, iso2c, iso3c, iso3n) %>%
  rename(country_id = vdem) 

# Find unique countries in the dataset
countries <- ds.basic %>% 
  filter(year ==2021) %>%
  dplyr::select(country, continent) %>%
  distinct()

# Create a dataframe with the additional years
new_years <- expand.grid(country = unique(countries$country), year = 2022:2023)

new_years <- new_years %>%
  mutate(continent = NA, 
         cown = NA, 
         country_id = NA, 
         iso2c = NA, 
         iso3c = NA, 
         iso3n = NA)

# Combine the new years with the existing dataset

ds.basic <- bind_rows(ds.basic, new_years) %>%
  arrange(country, year) %>%
  group_by(country) %>%
  fill(continent, cown, country_id, iso2c, iso3c, iso3n)

##################################################
#### Merge relevant V-Dem data ####
##################################################


decay_sum <- function(tm, vl, decay) {
  last_time <- 0
  current_sum <- 0
  sums <- numeric(length(vl))
  ldecay <- (1-decay)
  for (i in 1:length(vl)) {
    delta <- as.numeric(tm[i] - last_time)
    current_sum <- current_sum * (ldecay * delta) + vl[i]
    last_time <- tm[i]
    sums[i] <- current_sum
  }
  sums
}

annual_depricition_function  <- function(year, var, decay) {
  last_time <- first(year)
  current_sum <- 0
  sums <- numeric(length(var))
  ldecay <- (1-decay)
  for (i in 1:length(var)) {
    delta <- as.numeric(year[i] - last_time)
    current_sum <- ldecay * (current_sum + decay * var[i])
    last_time <- year[i]
    sums[i] <- current_sum
  }
  sums
}


## Functions ##

na_interpolation2 <- function(x, option = "linear", maxgap) {
  library(imputeTS)
  library(dplyr)
  
  total_not_missing <- sum(!is.na(x))
  
  # check there is sufficient data for na_interpolation 
  if(total_not_missing < 2) {x} 
  
  else
    
    # replace takes an input vector, a T/F vector & replacement value
  {replace(
    # input vector is interpolated data
    # this will impute leading/lagging NAs which we don't want 
    imputeTS::na_interpolation(x, option = option, maxgap = maxgap), 
    
    # create T/F vector for NAs,  
    is.na(na.approx(x, na.rm = FALSE)), 
    
    # replace TRUE with NA in input vector  
    NA) 
  }
}


vdem_subset_v14 <- vdem_full_v14 %>%
  group_by(country_id) %>%
  mutate(v2x_polyarchy_interpolated = na_interpolation2(v2x_polyarchy, option = "linear", maxgap = 5), 
         v2x_polyarchy_sd_interpolated = na_interpolation2(v2x_polyarchy_sd, option = "linear", maxgap = 5)) %>%
  drop_na(v2x_polyarchy_interpolated) %>%
  mutate(democracy_stock = annual_depricition_function(year, v2x_polyarchy_interpolated, 0.025), 
         democracy_stock_sd = annual_depricition_function(year, v2x_polyarchy_sd_interpolated, 0.025)) %>% # 2.5% annual depreciation rate
  dplyr::select(country_id, country_text_id, historical_date, year, v2x_polyarchy, v2x_polyarchy_sd, democracy_stock, democracy_stock_sd, e_democ, e_wb_pop)

summary(vdem_subset_v14$democracy_stock)
summary(vdem_subset_v14$democracy_stock_sd)

ds.basic <- ds.basic %>%
  left_join(vdem_subset_v14, by = c("country_id", "year"))

##############################################################
#### Merge relevant executive constraints data from V-Dem ####
##############################################################

vdem_subset_v14 <- vdem_full_v14 %>%
  group_by(country_id) %>%
  dplyr::select(country_id, year, starts_with("v2xlg_legcon"), starts_with("v2x_jucon"), starts_with("v2x_horacc"), 
                starts_with("v2juhcind"), starts_with("v2juncind"), starts_with("v2juhccomp"), 
                starts_with("v2jucomp"), starts_with("v2exrescon"), starts_with("v2lgotovst"), starts_with("v2lginvstp"),
                starts_with("v2lgbicam"), starts_with("v2lgqstexp"))

ds.basic <- ds.basic %>%
  left_join(vdem_subset_v14, by = c("country_id", "year"))

##################################################
#### Merge relevant V-Party data ####
##################################################

summary(v_party_v2$v2xpa_antiplural)


ggplot(v_party_v2, aes(x = v2xpa_antiplural)) +
 geom_density() +
 theme_clean()

# Construct standard Deviation for the Anti-pluralist Party Index, not SD available in original dataset but 68% uncertainty intervals
v_party_v2 <- v_party_v2 %>%
  mutate(v2xpa_antiplural_sd = (v2xpa_antiplural_codehigh-v2xpa_antiplural_codelow)/2) # 68% uncertainty intervals used by V-Dem 

# Comments: A 68% credible interval corresponds to approximately one standard deviation from the mean in both directions, 
# assuming the posterior distribution is approximately normal. Therefore, I used the width of the interval to estimate the 
# standard deviation. The sd was estimated by using the width of the uncertainty interval and dividing it by 1.96 

summary(v_party_v2$v2xpa_antiplural)
summary(v_party_v2$v2xpa_antiplural_sd)

## Tets if the assumtion works? # 
v_party_v2 <- v_party_v2 %>%
  mutate(v2xpa_antiplural_test_low = v2xpa_antiplural - v2xpa_antiplural_sd, 
         v2xpa_antiplural_test_high = v2xpa_antiplural + v2xpa_antiplural_sd)

test <- v_party_v2 %>%
  dplyr::select(country_text_id, year,v2xpa_antiplural, v2xpa_antiplural_sd,v2xpa_antiplural_codehigh, v2xpa_antiplural_test_high, 
                v2xpa_antiplural_codelow, v2xpa_antiplural_test_low) %>%
  mutate(diff_code_high = v2xpa_antiplural_codehigh-v2xpa_antiplural_test_high, 
         diff_code_low = v2xpa_antiplural_codelow - v2xpa_antiplural_test_low)

summary(test)
# the mean for the difference between the original 68% uncertainty intervals and the newly constructed confidence 
# intervals is near to 0. # the assumption of the normal distribution seems to hold

## Building country-party-year nested dataset from scratch ##

party_year_dataset <- countrycode::codelist_panel %>% 
  filter(year >= 1970) %>% # reduce to post-1946 period
  rename(country = country.name.en) %>% # reduce to post-1946 period
  dplyr::select(country, year, vdem) %>%
  rename(country_id = vdem)


## merging data from V-Party to country-party-year dataset ##

fill_na <- function(vector, n){
  if (n == 0) {
    vector
  } else {
    fill_na(
      vector = dplyr::coalesce(vector, dplyr::lag(vector)),
      n = n - 1
    )
  }
}


v_party_subset <- v_party_v2 %>%
  dplyr::select(country_id, year, v2paid, v2xpa_antiplural, v2xpa_antiplural_sd, v2paseatshare,
                v2pavote, v2paallian, v2pavallian, v2pagovsup) 

table(v_party_subset$v2pagovsup)

v_party_subset <- v_party_subset %>%
  mutate(gover_supp = ifelse(v2pagovsup %in% c(0,1,2), 1, 
                             ifelse(v2pagovsup %in% c(3), 0, NA)))

table(v_party_subset$gover_supp)

v_party_subset_country <- v_party_subset %>%
  mutate(v2paseatshare = v2paseatshare/100) %>%
  group_by(country_id, year) %>%
  summarize(anti_pluralist_party_index_gov = sum(v2xpa_antiplural[gover_supp==1]*v2paseatshare[gover_supp==1],na.rm = T), 
            anti_pluralist_party_index_gov_sd = sum(v2xpa_antiplural_sd[gover_supp==1]*v2paseatshare[gover_supp==1],na.rm = T), 
            anti_pluralist_party_index_opp = sum(v2xpa_antiplural[gover_supp==0]*v2paseatshare[gover_supp==0],na.rm = T), 
            anti_pluralist_party_index_opp_sd = sum(v2xpa_antiplural_sd[gover_supp==0]*v2paseatshare[gover_supp==0],na.rm = T)) %>%
  mutate(anti_pluralist_party_index = anti_pluralist_party_index_gov + anti_pluralist_party_index_opp, 
         anti_pluralist_party_index_sd = anti_pluralist_party_index_gov_sd + anti_pluralist_party_index_opp_sd, 
         anti_pluralist_party_index = ifelse(anti_pluralist_party_index ==0, NA, anti_pluralist_party_index), 
         anti_pluralist_party_index_sd = ifelse(anti_pluralist_party_index_sd == 0, NA, anti_pluralist_party_index_sd))
  
party_year_dataset <- party_year_dataset %>%
  left_join(v_party_subset_country, by = c("country_id", "year")) %>%
  fill(country_id) %>%
  mutate(anti_pluralist_party_index = fill_na(anti_pluralist_party_index, 10), 
         anti_pluralist_party_index_sd = fill_na(anti_pluralist_party_index_sd, 10))

ggplot(party_year_dataset, aes(x = anti_pluralist_party_index)) +
  geom_density() +
  theme_clean()

summary(party_year_dataset$anti_pluralist_party_index)
summary(party_year_dataset$anti_pluralist_party_index_sd)

party_year_dataset <- party_year_dataset %>%
  dplyr::select(-c(country_id, starts_with("anti_pluralist_party_index_gov"), starts_with("anti_pluralist_party_index_opp"))) %>%
  mutate(anti_pluralist_party_index_rev = scales::rescale(anti_pluralist_party_index,  to = c(1, 0.0098)))

summary(party_year_dataset$anti_pluralist_party_index_rev)

ggplot(party_year_dataset, aes(x = anti_pluralist_party_index_rev)) +
  geom_density() +
  theme_clean()


## merge anti-pluralist vote share data into ds.basic ##

ds.basic <- ds.basic %>%
  left_join(party_year_dataset, by = c("country", "year")) 


##################################################
#### Polarization                             ####
##################################################

vdem_subset_v14 <- vdem_full_v14 %>%
  dplyr::select(country_id, year, v2cacamps, v2cacamps_sd, v2cacamps_osp, v2cacamps_osp_sd, 
                v2caviol, v2caviol_sd, v2caviol_osp, v2caviol_osp_sd) 

summary(vdem_subset_v14$v2cacamps)
summary(vdem_subset_v14$v2cacamps_osp)
summary(vdem_subset_v14$v2caviol)
summary(vdem_subset_v14$v2caviol_osp)

vdem_subset_v14 <- vdem_subset_v14 %>%
  mutate(v2cacamps_rev = scales::rescale(v2cacamps, to = c(3.934, -3.806)), 
         v2cacamps_osp_rev = scales::rescale(v2cacamps_osp, to = c(3.994, 0.044)), 
         v2caviol_rev = scales::rescale(v2caviol, to = c(4.164, -3.461)), 
         v2caviol_osp_rev = scales::rescale(v2caviol_osp, to = c(3.99, 0.038)))

ds.basic <- ds.basic %>%
  left_join(vdem_subset_v14, by = c("country_id", "year"))

ggplot(vdem_subset_v14, aes(x = v2cacamps_rev)) +
  geom_density() +
  theme_clean()
ggplot(vdem_subset_v14, aes(x = v2cacamps_osp_rev)) +
  geom_density() +
  theme_clean()


##################################################
#### Robustness of Civil Society: V-Dem       ####
##################################################

vdem_subset_v14 <- vdem_full_v14 %>%
  dplyr::select(country_id, year, starts_with(c("v2xcs_ccsi", "v2cseeorgs", "v2csreprss", "v2csprtcpt")))

summary(vdem_subset_v14$v2xcs_ccsi)

ds.basic <- ds.basic %>%
  left_join(vdem_subset_v14, by = c("country_id", "year"))

vdem_subset_v14 <- vdem_full_v14 %>%
  dplyr::select(country_id, year, starts_with(c("v2x_cspart", "v2pscnslnl", "v2cscnsult", "v2csgender")))

summary(vdem_subset_v14$v2xcs_ccsi)

ds.basic <- ds.basic %>%
  left_join(vdem_subset_v14, by = c("country_id", "year"))


##################################################
#### Equal access index: V-Dem                ####
##################################################
#alternative to Vanhanen power resources measure that is available for years after 2000

vdem_subset_v14 <- vdem_full_v14 %>%
  dplyr::select(country_id, year, starts_with(c("v2xeg_eqaccess", "v2pepwrgen", "v2pepwrsoc", "v2pepwrses")))

summary(vdem_subset_v14$v2xeg_eqaccess)

ds.basic <- ds.basic %>%
  left_join(vdem_subset_v14, by = c("country_id", "year"))


##################################################
####Equal distribution of resources index: V-Dem                ####
##################################################
#alternative to Vanhanen power resources measure that is available for years after 2000

vdem_subset_v14 <- vdem_full_v14 %>%
  dplyr::select(country_id, year, starts_with(c("v2xeg_eqdr", "v2dlencmps", "v2dlunivl", "v2peedueq", "v2pehealth")))

summary(vdem_subset_v14$v2xeg_eqdr)

ds.basic <- ds.basic %>%
  left_join(vdem_subset_v14, by = c("country_id", "year"))


##################################################
#### Satisfaction with Democracy              ####
##################################################

democratic_satisfaction_data <- democratic_satisfaction_data %>%
  rename(iso3c = ISO_code, 
         year = Year) %>%
  dplyr::select(-Country)

ds.basic <- ds.basic %>%
  left_join(democratic_satisfaction_data, by = c("iso3c", "year"))


##################################################
#### Democratic Mood                          ####
##################################################

democratic_mood_data <- democratic_mood_data %>%
  rename(iso3c = ISO3c, 
         year = Year) %>%
  dplyr::select(-Country)

ds.basic <- ds.basic %>%
  left_join(democratic_mood_data, by = c("iso3c", "year"))


###################################################
#### Political Trust BJPS Study ###################
###################################################

civil_mood_data$iso3c <-  countrycode(civil_mood_data$Country, origin = "country.name", destination = "iso3c")
gov_mood_data$iso3c <-  countrycode(gov_mood_data$Country, origin = "country.name", destination = "iso3c")
leg_mood_data$iso3c <-  countrycode(leg_mood_data$Country, origin = "country.name", destination = "iso3c")
parl_mood_data$iso3c <-  countrycode(parl_mood_data$Country, origin = "country.name", destination = "iso3c")
police_mood_data$iso3c <-  countrycode(police_mood_data$Country, origin = "country.name", destination = "iso3c")
polpar_mood_data$iso3c <-  countrycode(polpar_mood_data$Country, origin = "country.name", destination = "iso3c")

civil_mood_data <- civil_mood_data %>% dplyr::select(-c(First_yr, Country))
gov_mood_data <- gov_mood_data %>% dplyr::select(-c(First_yr, Country))
leg_mood_data <- leg_mood_data %>% dplyr::select(-c(First_yr, Country))
parl_mood_data <- parl_mood_data %>% dplyr::select(-c(First_yr, Country))
police_mood_data <- police_mood_data %>% dplyr::select(-c(First_yr, Country))
polpar_mood_data <- polpar_mood_data %>% dplyr::select(-c(First_yr, Country))

ds.basic <- ds.basic %>%
  left_join(civil_mood_data, by = c("iso3c", "year" = "Year")) %>%
  left_join(gov_mood_data, by = c("iso3c", "year" = "Year")) %>%
  left_join(leg_mood_data, by = c("iso3c", "year" = "Year")) %>%
  left_join(parl_mood_data, by = c("iso3c", "year" = "Year")) %>%
  left_join(police_mood_data, by = c("iso3c", "year" = "Year")) %>%
  left_join(polpar_mood_data, by = c("iso3c", "year" = "Year")) 


summary(ds.basic$civil)

ds.basic <- ds.basic %>%
  group_by(country_text_id) %>%
  mutate(helper_civil = ifelse(sum(!is.na(civil)) > 1, TRUE, FALSE), 
         helper_gov = ifelse(sum(!is.na(gov)) > 1, TRUE, FALSE), 
         helper_leg = ifelse(sum(!is.na(leg)) > 1, TRUE, FALSE), 
         helper_parl = ifelse(sum(!is.na(parl)) > 1, TRUE, FALSE), 
         helper_police = ifelse(sum(!is.na(police)) > 1, TRUE, FALSE), 
         helper_polpar = ifelse(sum(!is.na(polpar)) > 1, TRUE, FALSE), 
         helper_satis = ifelse(sum(!is.na(Satis)) > 1, TRUE, FALSE)) %>%
  mutate(civil = ifelse(helper_civil == T, 
                             imputeTS::na_interpolation(civil, option = "spline", maxgap = 4),
                             civil), 
         civil_sd = ifelse(helper_civil == T, 
                        imputeTS::na_interpolation(civil_sd, option = "spline", maxgap = 4),
                        civil_sd), 
         gov = ifelse(helper_gov == T, 
                        imputeTS::na_interpolation(gov, option = "spline", maxgap = 4),
                      gov), 
         gov_sd = ifelse(helper_gov == T, 
                           imputeTS::na_interpolation(gov_sd, option = "spline", maxgap = 4),
                         gov_sd), 
         leg = ifelse(helper_leg == T, 
                        imputeTS::na_interpolation(leg, option = "spline", maxgap = 4),
                      leg), 
         leg_sd = ifelse(helper_leg == T, 
                           imputeTS::na_interpolation(leg_sd, option = "spline", maxgap = 4),
                         leg_sd), 
         parl = ifelse(helper_parl == T, 
                        imputeTS::na_interpolation(parl, option = "spline", maxgap = 4),
                       parl), 
         parl_sd = ifelse(helper_parl == T, 
                           imputeTS::na_interpolation(parl_sd, option = "spline", maxgap = 4),
                          parl_sd), 
         police = ifelse(helper_police == T, 
                        imputeTS::na_interpolation(police, option = "spline", maxgap = 4),
                        police), 
         police_sd = ifelse(helper_police == T, 
                           imputeTS::na_interpolation(police_sd, option = "spline", maxgap = 4),
                           police_sd), 
         polpar = ifelse(helper_polpar== T, 
                        imputeTS::na_interpolation(polpar, option = "spline", maxgap = 4),
                        polpar), 
         polpar_sd = ifelse(helper_polpar == T, 
                           imputeTS::na_interpolation(polpar_sd, option = "spline", maxgap = 4),
                           polpar_sd), 
         Satis = ifelse(helper_satis== T, 
                         imputeTS::na_interpolation(Satis, option = "spline", maxgap = 4),
                        Satis), 
         Satis_sd = ifelse(helper_satis == T, 
                            imputeTS::na_interpolation(Satis_sd, option = "spline", maxgap = 4),
                           Satis_sd))

################################################################################
################################################################################

#### Inspect Patterns (Missingness etc.) in the data ####

summary(ds.basic)

ds.basic.subset.2000 <- ds.basic %>%
  filter(year >=2000) %>%
  dplyr::select(-c(continent, cown, country_id, e_wb_pop))

prepare_missing_values_graph(ds.basic.subset.2000, ts_id = "year", no_factors = FALSE, binary = FALSE)

summary(ds.basic$v2cacamps)
summary(ds.basic$anti_pluralist_party_index)

ds.basic <- ds.basic %>%
  ungroup() %>%
  mutate(v2cacamps_rev = scales::rescale(v2cacamps, to = c(-3.934,3.806)), 
         anti_pluralist_party_index_rev = scales::rescale(anti_pluralist_party_index, to = c(1, 0)))
  
summary(ds.basic$v2cacamps_rev)
summary(ds.basic$anti_pluralist_party_index_rev)
summary(ds.basic$anti_pluralist_party_index)

saveRDS(ds.basic.subset.2000, "data/ds.basic.subset.2000.rds")
ds.basic.subset.2000 <- readRDS("data/ds.basic.subset.2000.rds")


################################################################################
################################################################################
# To estimate the different dimensions of resilience capacity, you need to run the R-script 01_BFA.R.
# This script approximately will run multiple hours using the rstan environment with following parameters:
# ITER <- 200
# BURNIN <- 10000
# MCMC <- 10000
# THIN <- 200

source("00_BFA.R")


################################################################################
################################################################################

# To run the Bayesian Factor Analysis to get the Resilience Capacity Components, please us the following script:
# 01_BFA.R
# This script will run the BFA analysis and store the posterior files as well as the thin_post files in the folder BFA. 
# What follows is the construction of the additive and multiplicative Resilience Capacity Index based on the posteriors from the BFA
# to include measurement error in the respective index. 

#####
# Code from Vutils from V-Dem

# load matrix
load_matrix <- function(m, drop.vignettes = TRUE) {
  stopifnot(is.matrix(m))
  if (isTRUE(drop.vignettes))
    m <- m[!is_vignette(rownames(m)),, drop = F]
  m
}

# bfa_stretch_z_sample
bfa_stretch_z_sample <- function(ll, utable) {
  stretch(load_matrix(ll[[1]]$thin_post), 
          ll[[1]]$country_dates, 
          rule_366 = TRUE, 
          interpolate = FALSE, 
          utable = utable)
}

hli_summary <- function(full.ma, VARNAME) {
  info("Summarising posteriors for " %^% VARNAME)
  cd <- dist_summary(t(full.ma))
  cy <- cy.day_mean(cd, historical_date, country_text_id)
  cd <- fix_stat_columns(cd, VARNAME)
  cy <- fix_stat_columns(cy, VARNAME)
  return(list(cd = cd, cy = cy, thin_post = full.ma))
}

hpd <- function(m, prob) {
  out <- coda::HPDinterval(coda::as.mcmc(m), prob = prob)
  
  suffix <- sub("^.*[.]", "", prob)
  colnames(out) <- c("codelow", "codehigh") %^% suffix
  
  out
}

colMedians <- function(x) UseMethod("colMedians")
colSDs <- function(x) UseMethod("colSDs")
rowMedians <- function(x) UseMethod("rowMedians")
rowSDs <- function(x) UseMethod("rowSDs")

b_summary <- function(m) {
  
  m_mean <- colMeans(m)
  m_median <- colMedians(m)
  m_sd <- colSDs(m)
  hpd95 <- hpd(m, .95)
  hpd68 <- hpd(m, .68)
  
  out <- cbind(m_mean, m_median, m_sd, hpd68, hpd95, colnames(m))
  colnames(out) <- c("mean", "median", "sd", "codelow68", "codehigh68",
                     "codelow95", "codehigh95", "coder_id")
  rownames(out) <- NULL
  out
}


##########

## Load Posterior files 

UTABLE <- "../data/bfa/country_unit.rds"
INDIR  <- "../data/bfa/results/thin_post/"
OUTDIR <- "../data/bfa/results"
ITER <- 900

VARS <- c("thin_post_institutional_dimension", "thin_post_political_parties_dimension", "thin_post_civil_society_dimension", 
          "thin_post_political_community_dimension")

utable <- readRDS(UTABLE)
utable_names <- with(utable, paste(country_text_id, historical_date))

# To preserve gaps we need to find where there's jumps in our
# utable. Extract those dates to insert as NA in our individual C vars
# so that when we stretch we don't fill across those periods.
ll <- with(utable, split(year, country_text_id)) %>%
  lapply(function(v) {
    v <- sort(v)
    gapstart <- v[c(diff(v) > 1, F)]
    
    if (length(gapstart) > 0)
      as.Date((gapstart + 1) %^% "-12-31")
  }) %>% Filter(Negate(is.null), .)

gap_dates <- lapply(names(ll), function(s) paste(s, ll[[s]])) %>% unlist

# Load each matrix into a list
vars.ll <- lapply(setNames(VARS, VARS), function(varname) {
  # Construct the file path
  file_path <- file.path(INDIR, paste0(varname, ".rds"))
  
  # Read the matrix
  m <- readRDS(file_path)
  
  # Set column names
  colnames(m) <- rep(varname, ncol(m))
  
  # Remove vignette rows and select columns based on ITER
  # This assumes you want to subset into equally sized groups defined by ITER
  out <- m[!is_vignette(rownames(m)), seq(1, by = ncol(m) / ITER, length.out = ITER)]
  
  return(out)
})

###
# Start by conforming the dimensions b/w all matrices.
cvar_names <- Reduce(union, lapply(vars.ll, rownames))
vars.ll <- lapply(vars.ll, vutils::stretch, cvar_names, gaps = F, preserve.na = T)

# We'll need this after running the model to set the rows we want
# missing (50% of input vars were missing) to NA so that we don't
# stretch over those dates
missing.ma <- lapply(vars.ll, function(ma) rowSums(is.na(ma)) == ncol(ma)) %>%
  do.call(rbind, .)

dropped_dates <- colnames(missing.ma)[colMeans(missing.ma) > .5]

if (length(dropped_dates) > 0) {
  sprintf("Found %d country-dates with >50%% missingness",
          length(dropped_dates)) %>% info
}

vars.ll <- lapply(vars.ll, function(m) m[!rownames(m) %in% dropped_dates,, drop = F])
cvar_names <- setdiff(cvar_names, dropped_dates)

stopifnot(do.call(all_identical, lapply(vars.ll, rownames)))

###
# Initial values. n_obs = total output length of `xi`, i.e. total
# number of unique country-dates.
n_pars <- length(VARS)
n_obs <- length(cvar_names)

sprintf("Found %d total obs", n_obs) %>% info

# ==========================================================================
# Generate raw posterior files for additive ResCap, multiplicative ResCap and ResCap.
# ==========================================================================

# Functions
# --------------------------------------------------------------------------
calc_ResCap <- function(thin_post_institutional_dimension, thin_post_political_parties_dimension, 
                        thin_post_civil_society_dimension, thin_post_political_community_dimension) {
  
  # multiplicative ResCap
  # -- every factor multiplies
  post.mrescap <-
    thin_post_institutional_dimension * thin_post_political_parties_dimension * thin_post_civil_society_dimension * thin_post_political_community_dimension
 
   # additive ResCap
  # -- every factor adds
  post.arescap <-
    0.25 * thin_post_institutional_dimension + 0.25 * thin_post_political_parties_dimension + 0.25 * thin_post_civil_society_dimension +
    0.25 * thin_post_political_community_dimension 
  
  # (Combined) ResCap Index
  post.rescap <- 0.5 * post.mrescap + 0.5 * post.arescap
  
  list(post.mrescap, post.arescap, post.rescap)
  
}

# Apply calc_ResCap to each column
output <- lapply(1:ncol(vars.ll[[1]]), function(idx) {
  calc_ResCap(
    vars.ll$thin_post_institutional_dimension[, idx, drop = FALSE],
    vars.ll$thin_post_political_parties_dimension[, idx, drop = FALSE],
    vars.ll$thin_post_civil_society_dimension[, idx, drop = FALSE],
    vars.ll$thin_post_political_community_dimension[, idx, drop = FALSE]
  )
})

# now we have a output list with 900 elements for each of the three ResCap indices (MResCap, AResCap, ResCap)
# we thus need to sum up this 900 elements from the posteriors to to point estimates, sd, and uncertainty intervals

rownames_output <- rownames(rownames_output <- output[[1]][[1]])

matrix_mrescap <- sapply(output, function(sublist) {
  # Extract the desired element (assuming it's convertible to numeric or intended type)
  sublist[[1]]
})

rownames(matrix_mrescap) <- rownames_output

matrix_arescap <- sapply(output, function(sublist) {
  # Extract the desired element (assuming it's convertible to numeric or intended type)
  sublist[[2]]
})
rownames(matrix_arescap) <- rownames_output


matrix_rescap <- sapply(output, function(sublist) {
  # Extract the desired element (assuming it's convertible to numeric or intended type)
  sublist[[3]]
})
rownames(matrix_rescap) <- rownames_output

utable_names <- with(utable, paste(country_text_id, historical_date))
full_names <- intersect(utable_names, rownames_output)

#################################################################################
##################################################################################
#one can use dist_summary() with t() function (transpose) to generate the final country-year datasets 

mrescap_final_cy.df <- dist_summary(t(matrix_mrescap), full_names)
arescap_final_cy.df <- dist_summary(t(matrix_arescap), full_names)
rescap_final_cy.df <- dist_summary(t(matrix_rescap), full_names)

mrescap_final_cy.df <- mrescap_final_cy.df %>%
  rename(MResCap = median, 
         MResCap_sd = sd, 
         MResCap_codelow = codelow68, 
         MResCap_codehigh = codehigh68) %>%
  dplyr::select(-c(mean, codelow95, codehigh95)) %>%
  mutate(year = year(historical_date)) %>%
  dplyr::select(country_text_id, year, everything(), -historical_date) 

arescap_final_cy.df <- arescap_final_cy.df %>%
  rename(AResCap = median, 
         AResCap_sd = sd, 
         AResCap_codelow = codelow68, 
         AResCap_codehigh = codehigh68) %>%
  dplyr::select(-c(mean, codelow95, codehigh95)) %>%
  mutate(year = year(historical_date)) %>%
  dplyr::select(country_text_id, year, everything(), -historical_date) 

rescap_final_cy.df <- rescap_final_cy.df %>%
  rename(ResCap = median, 
         ResCap_sd = sd, 
         ResCap_codelow = codelow68, 
         ResCap_codehigh = codehigh68) %>%
  dplyr::select(-c(mean, codelow95, codehigh95)) %>%
  mutate(year = year(historical_date)) %>%
  dplyr::select(country_text_id, year, everything(), -historical_date) 

ds.basic.subset.2000 <- ds.basic.subset.2000 %>%
  left_join(rescap_final_cy.df, by = c("country_text_id", "year")) %>%
  left_join(arescap_final_cy.df, by = c("country_text_id", "year")) %>%
  left_join(mrescap_final_cy.df, by = c("country_text_id", "year"))

#### load four dimensions and there point estimates as well as uncertainty intervals 

institutional_dimension <- readRDS("../data/bfa/results/results/cy/bayes.institutional_dimension_cy.rds")
political_parties_dimension <- readRDS("../data/bfa/results/results/cy/bayes.political_parties_dimension_cy.rds")
civil_society_dimension<- readRDS("../data/bfa/results/results/cy/bayes.civil_society_dimension_cy.rds")
political_community_dimension <- readRDS("../data/bfa/results/results/cy/bayes.political_community_dimension_cy.rds")

institutional_dimension <- institutional_dimension %>%
  rename(inst_dim = median, 
         inst_dim_sd = sd, 
         inst_dim_codelow = codelow68, 
         inst_dim_codehigh = codehigh68) %>%
  dplyr::select(-c(mean, codelow95, codehigh95)) %>%
  dplyr::select(country_text_id, year, everything()) 

political_parties_dimension <- political_parties_dimension %>%
  rename(polpar_dim = median, 
         polpar_dim_sd = sd, 
         polpar_dim_codelow = codelow68, 
         polpar_dim_codehigh = codehigh68) %>%
  dplyr::select(-c(mean, codelow95, codehigh95)) %>%
  dplyr::select(country_text_id, year, everything()) 

civil_society_dimension <- civil_society_dimension %>%
  rename(civsoc_dim = median, 
         civsoc_dim_sd = sd, 
         civsoc_dim_codelow = codelow68, 
         civsoc_dim_codehigh = codehigh68) %>%
  dplyr::select(-c(mean, codelow95, codehigh95)) %>%
  dplyr::select(country_text_id, year, everything()) 


political_community_dimension <- political_community_dimension %>%
  rename(polcom_dim = median, 
         polcom_dim_sd = sd, 
         polcom_dim_codelow = codelow68, 
         polcom_dim_codehigh = codehigh68) %>%
  dplyr::select(-c(mean, codelow95, codehigh95)) %>%
  dplyr::select(country_text_id, year, everything()) 


ds.basic.subset.2000 <- ds.basic.subset.2000 %>%
  left_join(institutional_dimension, by = c("country_text_id", "year")) %>%
  left_join(political_parties_dimension, by = c("country_text_id", "year")) %>%
  left_join(civil_society_dimension, by = c("country_text_id", "year")) %>%
  left_join(political_community_dimension, by = c("country_text_id", "year")) 

################################################################################
################################################################################

#### Figure 2 Main Paper  ####

vdem_full_v14 <- readRDS(file.path("data/vdem14", "V-Dem-CY-Full+Others-v14.rds"))

vdem_subset_v14 <- vdem_full_v14 %>%
  dplyr::select(country_text_id, year, COWcode)

bmr_data <- read.csv("data/bmr/bmr_data_extended.csv")
bmr_data <- bmr_data %>% dplyr::select(-X)

bmr_data$country_text_id <- countrycode::countrycode(bmr_data$country, origin = "country.name", destination = "iso3c")

ds.basic.subset.2000_subset_democracies <- ds.basic.subset.2000 %>%
  left_join(vdem_subset_v14, by = c("country_text_id", "year")) %>%
  left_join(bmr_data, by = c("country_text_id", "year")) %>%
  filter(democracy == 1)

ds.basic.subset.2000_subset_democracies <- ds.basic.subset.2000_subset_democracies %>% distinct(country_id, year, .keep_all = TRUE)


f1 <- ggplot(ds.basic.subset.2000_subset_democracies, aes(x = ResCap)) +
  geom_density(fill="gray")+
  geom_vline(aes(xintercept=mean(ResCap, na.rm = TRUE)), color="blue",
             linetype="dashed")+
  labs(title="DRC Index",x="DRC Index", y = "Density")+
  theme_classic() +
  xlim(0,1)

f2 <- ggplot(ds.basic.subset.2000_subset_democracies, aes(x = AResCap)) +
  geom_density(fill="gray")+
  geom_vline(aes(xintercept=mean(AResCap, na.rm = TRUE)), color="blue",
             linetype="dashed")+
  labs(title="additive DRC Index",x="additive DRC Index", y = "Density")+
  theme_classic() +
  xlim(0,1)

f3 <- ggplot(ds.basic.subset.2000_subset_democracies, aes(x = MResCap)) +
  geom_density(fill="gray")+
  geom_vline(aes(xintercept=mean(MResCap, na.rm = TRUE)), color="blue",
             linetype="dashed")+
  labs(title="multiplicative DRC Index",x="multiplicative DRC Index", y = "Density")+
  theme_classic() +
  xlim(0,1)

f4 <-  ggplot(ds.basic.subset.2000_subset_democracies, aes(x = MResCap, y = AResCap)) +
  geom_point(aes(alpha = 0.05)) +
  labs(title="",x="multiplicative DRC Index", y = "additive DRC Index")+
  theme_classic() +
  guides(alpha = "none")


ggarrange(f1, f2, f3, f4,
          labels = c("A", "B", "C", "D"),
          ncol = 2, nrow = 2)

ggsave("..//Figure_2.png", width = 30,
       height = 30,
       units = c("cm"), dpi = 1000)

ggsave("..//Figure_2_low_res.png", width = 30,
       height = 30,
       units = c("cm"), dpi = 400)

#### Figure 3 Main Paper ####

GER_F <- ggplot(subset(ds.basic.subset.2000, country %in% c("Germany")), 
                aes(x=year)) +
  geom_line(aes(y = ResCap), color = "black") +
  geom_ribbon(aes(ymin = ResCap_codelow, ymax = ResCap_codehigh), fill = "black", alpha = 0.2) +
  geom_line(aes(y = AResCap), color = "#00AFBB") +
  geom_ribbon(aes(ymin = AResCap_codelow, ymax = AResCap_codehigh),  fill = "#00AFBB", alpha = 0.2) +
  geom_line(aes(y = MResCap), color = "#E7B800") +
  geom_ribbon(aes(ymin = MResCap_codelow, ymax = MResCap_codehigh), fill = "#E7B800" , alpha = 0.2) +
  theme_pubclean()+
  ylim(0,1) + ggtitle("Germany") + theme(plot.title = element_text(hjust = 0.5)) +ylab("DRC Index")

NOR_F <- ggplot(subset(ds.basic.subset.2000, country %in% c("Norway")), 
                aes(x=year)) +
  geom_line(aes(y = ResCap), color = "black") +
  geom_ribbon(aes(ymin = ResCap_codelow, ymax = ResCap_codehigh), fill = "black", alpha = 0.2) +
  geom_line(aes(y = AResCap), color = "#00AFBB") +
  geom_ribbon(aes(ymin = AResCap_codelow, ymax = AResCap_codehigh),  fill = "#00AFBB", alpha = 0.2) +
  geom_line(aes(y = MResCap), color = "#E7B800") +
  geom_ribbon(aes(ymin = MResCap_codelow, ymax = MResCap_codehigh), fill = "#E7B800" , alpha = 0.2) +
  theme_pubclean()+
  ylim(0,1) + ggtitle("Norway") + theme(plot.title = element_text(hjust = 0.5)) +ylab("DRC Index")

RUS_F <- ggplot(subset(ds.basic.subset.2000, country %in% c("Russia")), 
                aes(x=year)) +
  geom_line(aes(y = ResCap), color = "black") +
  geom_ribbon(aes(ymin = ResCap_codelow, ymax = ResCap_codehigh), fill = "black", alpha = 0.2) +
  geom_line(aes(y = AResCap), color = "#00AFBB") +
  geom_ribbon(aes(ymin = AResCap_codelow, ymax = AResCap_codehigh),  fill = "#00AFBB", alpha = 0.2) +
  geom_line(aes(y = MResCap), color = "#E7B800") +
  geom_ribbon(aes(ymin = MResCap_codelow, ymax = MResCap_codehigh), fill = "#E7B800" , alpha = 0.2) +
  theme_pubclean()+
  ylim(0,1) + ggtitle("Russia") + theme(plot.title = element_text(hjust = 0.5)) +ylab("DRC Index")

USA_F <- ggplot(subset(ds.basic.subset.2000, country %in% c("United States")), 
                aes(x=year)) +
  geom_line(aes(y = ResCap), color = "black") +
  geom_ribbon(aes(ymin = ResCap_codelow, ymax = ResCap_codehigh), fill = "black", alpha = 0.2) +
  geom_line(aes(y = AResCap), color = "#00AFBB") +
  geom_ribbon(aes(ymin = AResCap_codelow, ymax = AResCap_codehigh),  fill = "#00AFBB", alpha = 0.2) +
  geom_line(aes(y = MResCap), color = "#E7B800") +
  geom_ribbon(aes(ymin = MResCap_codelow, ymax = MResCap_codehigh), fill = "#E7B800" , alpha = 0.2) +
  theme_pubclean()+
  ylim(0,1) + ggtitle("United States") + theme(plot.title = element_text(hjust = 0.5)) +ylab("DRC Index")


FRA_F <- ggplot(subset(ds.basic.subset.2000, country %in% c("France")), 
                aes(x=year)) +
  geom_line(aes(y = ResCap), color = "black") +
  geom_ribbon(aes(ymin = ResCap_codelow, ymax = ResCap_codehigh), fill = "black", alpha = 0.2) +
  geom_line(aes(y = AResCap), color = "#00AFBB") +
  geom_ribbon(aes(ymin = AResCap_codelow, ymax = AResCap_codehigh),  fill = "#00AFBB", alpha = 0.2) +
  geom_line(aes(y = MResCap), color = "#E7B800") +
  geom_ribbon(aes(ymin = MResCap_codelow, ymax = MResCap_codehigh), fill = "#E7B800" , alpha = 0.2) +
  theme_pubclean()+
  ylim(0,1) + ggtitle("France") + theme(plot.title = element_text(hjust = 0.5)) +ylab("DRC Index")

SOK_F <- ggplot(subset(ds.basic.subset.2000, country %in% c("South Korea")), 
                aes(x=year)) +
  geom_line(aes(y = ResCap), color = "black") +
  geom_ribbon(aes(ymin = ResCap_codelow, ymax = ResCap_codehigh), fill = "black", alpha = 0.2) +
  geom_line(aes(y = AResCap), color = "#00AFBB") +
  geom_ribbon(aes(ymin = AResCap_codelow, ymax = AResCap_codehigh),  fill = "#00AFBB", alpha = 0.2) +
  geom_line(aes(y = MResCap), color = "#E7B800") +
  geom_ribbon(aes(ymin = MResCap_codelow, ymax = MResCap_codehigh), fill = "#E7B800" , alpha = 0.2) +
  theme_pubclean()+
  ylim(0,1) + ggtitle("South Korea") + theme(plot.title = element_text(hjust = 0.5)) +ylab("DRC Index")

(USA_F + RUS_F + SOK_F) / (FRA_F + GER_F + NOR_F)

ggsave("..//Figure_3.png", width = 30,
       height = 20,
       units = c("cm"), dpi = 1000)

ggsave("..//Figure_3_low_res.png", width = 30,
       height = 20,
       units = c("cm"), dpi = 400)

#### Figure A1 Appendix  ####

f1 <- ggplot(ds.basic.subset.2000, aes(x = ResCap)) +
  geom_density(fill="gray")+
  geom_vline(aes(xintercept=mean(ResCap, na.rm = TRUE)), color="blue",
             linetype="dashed")+
  labs(title="DRC Index",x="DRC Index", y = "Density")+
  theme_classic() +
  xlim(0,1)

f2 <- ggplot(ds.basic.subset.2000, aes(x = AResCap)) +
  geom_density(fill="gray")+
  geom_vline(aes(xintercept=mean(AResCap, na.rm = TRUE)), color="blue",
             linetype="dashed")+
  labs(title="additive DRC Index",x="additive DRC Index", y = "Density")+
  theme_classic() +
  xlim(0,1)

f3 <- ggplot(ds.basic.subset.2000, aes(x = MResCap)) +
  geom_density(fill="gray")+
  geom_vline(aes(xintercept=mean(MResCap, na.rm = TRUE)), color="blue",
             linetype="dashed")+
  labs(title="multiplicative DRC Index",x="multiplicative DRC Index", y = "Density")+
  theme_classic() +
  xlim(0,1)

f4 <-  ggplot(ds.basic.subset.2000, aes(x = MResCap, y = AResCap)) +
  geom_point(aes(alpha = 0.05)) +
  labs(title="",x="multiplicative DRC Index", y = "additive DRC Index")+
  theme_classic() +
  guides(alpha = "none")


ggarrange(f1, f2, f3, f4,
          labels = c("A", "B", "C", "D"),
          ncol = 2, nrow = 2)

ggsave("..//Figure_A1.png", width = 30,
       height = 30,
       units = c("cm"), dpi = 1000)

ggsave("..//Figure_A1_low_res.png", width = 30,
       height = 30,
       units = c("cm"), dpi = 400)

#### World Average for all Democracies?

world_averages <- ds.basic.subset.2000 %>%
  left_join(bmr_data, by = c("iso3c" ="abbreviation_undp", "year")) %>%
  filter(democracy == 1) %>%
  group_by(year) %>%
  summarize(mean = mean(ResCap, na.rm = T), 
            mean_codelow = mean(ResCap_codelow, na.rm = T), 
            mean_codehigh = mean(ResCap_codehigh, na.rm = T))
  
world_averages_plot <- ggplot(world_averages, aes(x=year)) +
  geom_line(aes(y = mean), color = "#00AFBB") +
  geom_ribbon(aes(ymin = mean_codelow, ymax = mean_codehigh), fill = "#00AFBB", alpha = 0.2) +
  theme_pubclean()+
  ylim(0,1) + ggtitle("Wold Averages") + theme(plot.title = element_text(hjust = 0.5))

ggsave(plot = world_averages_plot, "..//world_averages.png", width = 30,
       height = 20,
       units = c("cm"), dpi = 1000)

ggsave(plot = world_averages_plot, "..//world_averages_low_res.png", width = 30,
       height = 20,
       units = c("cm"), dpi = 400)

#### Figure 4: Dimensions of Resilience Capacity ####

GER_F <- ggplot(subset(ds.basic.subset.2000, country %in% c("Germany")), 
                aes(x=year)) +
  geom_line(aes(y = inst_dim), color = "#A62631FF") +
  geom_ribbon(aes(ymin = inst_dim_codelow, ymax = inst_dim_codehigh), fill = "#A62631FF", alpha = 0.2) +
  geom_line(aes(y = polpar_dim), color = "#00AFBB") +
  geom_ribbon(aes(ymin = polpar_dim_codelow, ymax = polpar_dim_codehigh),  fill = "#00AFBB", alpha = 0.2) +
  geom_line(aes(y = civsoc_dim), color = "#E7B800") +
  geom_ribbon(aes(ymin = civsoc_dim_codelow, ymax = civsoc_dim_codehigh), fill = "#E7B800" , alpha = 0.2) +
  geom_line(aes(y = polcom_dim), color = "#272A59FF") +
  geom_ribbon(aes(ymin = polcom_dim_codelow, ymax = polcom_dim_codehigh), fill = "#272A59FF" , alpha = 0.2) +
  theme_pubclean()+
  ylim(0,1) + ggtitle("Germany") + theme(plot.title = element_text(hjust = 0.5)) + ylab("")

NOR_F <- ggplot(subset(ds.basic.subset.2000, country %in% c("Norway")), 
                aes(x=year)) +
  geom_line(aes(y = inst_dim), color = "#A62631FF") +
  geom_ribbon(aes(ymin = inst_dim_codelow, ymax = inst_dim_codehigh), fill = "#A62631FF", alpha = 0.2) +
  geom_line(aes(y = polpar_dim), color = "#00AFBB") +
  geom_ribbon(aes(ymin = polpar_dim_codelow, ymax = polpar_dim_codehigh),  fill = "#00AFBB", alpha = 0.2) +
  geom_line(aes(y = civsoc_dim), color = "#E7B800") +
  geom_ribbon(aes(ymin = civsoc_dim_codelow, ymax = civsoc_dim_codehigh), fill = "#E7B800" , alpha = 0.2) +
  geom_line(aes(y = polcom_dim), color = "#272A59FF") +
  geom_ribbon(aes(ymin = polcom_dim_codelow, ymax = polcom_dim_codehigh), fill = "#272A59FF" , alpha = 0.2) +
  theme_pubclean()+
  ylim(0,1) + ggtitle("Norway") + theme(plot.title = element_text(hjust = 0.5)) + ylab("")

RUS_F <- ggplot(subset(ds.basic.subset.2000, country %in% c("Russia")), 
                aes(x=year)) +
  geom_line(aes(y = inst_dim), color = "#A62631FF") +
  geom_ribbon(aes(ymin = inst_dim_codelow, ymax = inst_dim_codehigh), fill = "#A62631FF", alpha = 0.2) +
  geom_line(aes(y = polpar_dim), color = "#00AFBB") +
  geom_ribbon(aes(ymin = polpar_dim_codelow, ymax = polpar_dim_codehigh),  fill = "#00AFBB", alpha = 0.2) +
  geom_line(aes(y = civsoc_dim), color = "#E7B800") +
  geom_ribbon(aes(ymin = civsoc_dim_codelow, ymax = civsoc_dim_codehigh), fill = "#E7B800" , alpha = 0.2) +
  geom_line(aes(y = polcom_dim), color = "#272A59FF") +
  geom_ribbon(aes(ymin = polcom_dim_codelow, ymax = polcom_dim_codehigh), fill = "#272A59FF" , alpha = 0.2) +
  theme_pubclean()+
  ylim(0,1) + ggtitle("Russia") + theme(plot.title = element_text(hjust = 0.5)) + ylab("")
 
USA_F <- ggplot(subset(ds.basic.subset.2000, country %in% c("United States")), 
                aes(x=year)) +
  geom_line(aes(y = inst_dim), color = "#A62631FF") +
  geom_ribbon(aes(ymin = inst_dim_codelow, ymax = inst_dim_codehigh), fill = "#A62631FF", alpha = 0.2) +
  geom_line(aes(y = polpar_dim), color = "#00AFBB") +
  geom_ribbon(aes(ymin = polpar_dim_codelow, ymax = polpar_dim_codehigh),  fill = "#00AFBB", alpha = 0.2) +
  geom_line(aes(y = civsoc_dim), color = "#E7B800") +
  geom_ribbon(aes(ymin = civsoc_dim_codelow, ymax = civsoc_dim_codehigh), fill = "#E7B800" , alpha = 0.2) +
  geom_line(aes(y = polcom_dim), color = "#272A59FF") +
  geom_ribbon(aes(ymin = polcom_dim_codelow, ymax = polcom_dim_codehigh), fill = "#272A59FF" , alpha = 0.2) +
  theme_pubclean()+
  ylim(0,1) + ggtitle("United States") + theme(plot.title = element_text(hjust = 0.5)) + ylab("")


FRA_F <- ggplot(subset(ds.basic.subset.2000, country %in% c("France")), 
                aes(x=year)) +
  geom_line(aes(y = inst_dim), color = "#A62631FF") +
  geom_ribbon(aes(ymin = inst_dim_codelow, ymax = inst_dim_codehigh), fill = "#A62631FF", alpha = 0.2) +
  geom_line(aes(y = polpar_dim), color = "#00AFBB") +
  geom_ribbon(aes(ymin = polpar_dim_codelow, ymax = polpar_dim_codehigh),  fill = "#00AFBB", alpha = 0.2) +
  geom_line(aes(y = civsoc_dim), color = "#E7B800") +
  geom_ribbon(aes(ymin = civsoc_dim_codelow, ymax = civsoc_dim_codehigh), fill = "#E7B800" , alpha = 0.2) +
  geom_line(aes(y = polcom_dim), color = "#272A59FF") +
  geom_ribbon(aes(ymin = polcom_dim_codelow, ymax = polcom_dim_codehigh), fill = "#272A59FF" , alpha = 0.2) +
  theme_pubclean()+
  ylim(0,1) + ggtitle("France") + theme(plot.title = element_text(hjust = 0.5)) + ylab("")

SOK_F <- ggplot(subset(ds.basic.subset.2000, country %in% c("South Korea")), 
                aes(x=year)) +
  geom_line(aes(y = inst_dim), color = "#A62631FF") +
  geom_ribbon(aes(ymin = inst_dim_codelow, ymax = inst_dim_codehigh), fill = "#A62631FF", alpha = 0.2) +
  geom_line(aes(y = polpar_dim), color = "#00AFBB") +
  geom_ribbon(aes(ymin = polpar_dim_codelow, ymax = polpar_dim_codehigh),  fill = "#00AFBB", alpha = 0.2) +
  geom_line(aes(y = civsoc_dim), color = "#E7B800") +
  geom_ribbon(aes(ymin = civsoc_dim_codelow, ymax = civsoc_dim_codehigh), fill = "#E7B800" , alpha = 0.2) +
  geom_line(aes(y = polcom_dim), color = "#272A59FF") +
  geom_ribbon(aes(ymin = polcom_dim_codelow, ymax = polcom_dim_codehigh), fill = "#272A59FF" , alpha = 0.2) +
  theme_pubclean()+
  ylim(0,1) + ggtitle("South Korea") + theme(plot.title = element_text(hjust = 0.5)) + ylab("")

(USA_F + RUS_F + SOK_F) / (FRA_F + GER_F + NOR_F)

ggsave("..//Figure_4.png", width = 30,
       height = 20,
       units = c("cm"), dpi = 1000)

ggsave("..//Figure_4_low_res.png", width = 30,
       height = 20,
       units = c("cm"), dpi = 400)

ds.basic.subset.2000_subset <- ds.basic.subset.2000 %>%
  dplyr::select(country_text_id, country, country_id, year, iso2c, iso3c, iso3n, historical_date, 
                starts_with(c("ResCap", "MResCap", "AResCap", "inst_dim", "polpar_dim", 
                              "civsoc_dim", "polcom_dim")))
saveRDS(ds.basic.subset.2000_subset, "data/ds.basic.subset.2000_subset.RDS")
write_csv(ds.basic.subset.2000_subset, "data/ds.basic.subset.2000_subset.csv")

ds.basic.subset.2000_subset <- readRDS("data/ds.basic.subset.2000_subset.RDS")

writexl::write_xlsx(ds.basic.subset.2000_subset, "data/ds.basic.subset.2000_subset.xlsx")

################################################################################
################################################################################

#### Face Validity Test ####

figure2_data <- ds.basic.subset.2000 %>%
  dplyr::select(country, year, ResCap, ResCap_codelow, ResCap_codehigh) %>% 
  group_by(country) %>%
  mutate(helper_ResCap = ifelse(sum(is.na(ResCap)) >10, TRUE, FALSE)) %>%
  ungroup() %>%
  filter(helper_ResCap == FALSE) %>%
  dplyr::select(-helper_ResCap) 

figure2_data_country <- unique(figure2_data$country)

figure2_data_country_A <- figure2_data_country[1:30]
figure2_data_country_B <- figure2_data_country[31:60]
figure2_data_country_C <- figure2_data_country[61:90]
figure2_data_country_D <- figure2_data_country[91:127]


## Figure A2 ##
ggplot(subset(figure2_data, country %in% figure2_data_country_A), 
       aes(x=year, y= ResCap, ymin = ResCap_codelow, ymax = ResCap_codehigh)) +
  geom_line(color = "#E7B800") +
  geom_ribbon(fill ="#E7B800", alpha = 0.3) +
  facet_wrap(~country) +
  theme_pubclean() +
  scale_x_continuous(name="Year", breaks=seq(2000,2023,10)) +
  ylim(0,1)+ ylab("DRC Index")

ggsave("..//Figure_A2.png", width = 20,
       height = 30,
       units = c("cm"), dpi = 400)

## Figure A3 ##
ggplot(subset(figure2_data, country %in% figure2_data_country_B), 
       aes(x=year, y= ResCap, ymin = ResCap_codelow, ymax = ResCap_codehigh)) +
  geom_line(color = "#E7B800") +
  geom_ribbon(fill ="#E7B800", alpha = 0.3) +
  facet_wrap(~country) +
  theme_pubclean() +
  scale_x_continuous(name="Year", breaks=seq(2000,2023,10)) +
  ylim(0,1)+ ylab("DRC Index")

ggsave("..//Figure_A3.png", width = 20,
       height = 30,
       units = c("cm"), dpi = 400)

## Figure A4 ##
ggplot(subset(figure2_data, country %in% figure2_data_country_C), 
       aes(x=year, y= ResCap, ymin = ResCap_codelow, ymax = ResCap_codehigh)) +
  geom_line(color = "#E7B800") +
  geom_ribbon(fill ="#E7B800", alpha = 0.3) +
  facet_wrap(~country) +
  theme_pubclean() +
  scale_x_continuous(name="Year", breaks=seq(2000,2023,10)) +
  ylim(0,1)+ ylab("DRC Index")

ggsave("..//Figure_A4.png", width = 20,
       height = 30,
       units = c("cm"), dpi = 400)

## Figure A5 ##
ggplot(subset(figure2_data, country %in% figure2_data_country_D), 
       aes(x=year, y= ResCap, ymin = ResCap_codelow, ymax = ResCap_codehigh)) +
  geom_line(color = "#E7B800") +
  geom_ribbon(fill ="#E7B800", alpha = 0.3) +
  facet_wrap(~country) +
  theme_pubclean() +
  scale_x_continuous(name="Year", breaks=seq(2000,2023,10)) +
  ylim(0,1) + ylab("DRC Index")

ggsave("../results/Figure_A5.png", width = 20,
       height = 30,
       units = c("cm"), dpi = 400)